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We consider the incompressible Navier-Stokes equations with periodic boundary con¬ 
ditions and time-independent forcing. For this type of flow, we derive adjoint equations 
whose trajectories converge asymptotically to the equilibrium and traveling wave solu¬ 
tions of the Navier-Stokes equations. Using the adjoint equations, arbitrary initial condi¬ 
tions evolve to the vicinity of a (relative) equilibrium at which point a few Newton-type 
iterations yield the desired (relative) equilibrium solution. We apply this adjoint-based 
method to a chaotic two-dimensional Kolmogorov flow. A convergence rate of 100% is 
observed, leading to the discovery of 21 new steady state and traveling wave solutions at 
Reynolds number Re = 40. Some of the new invariant solutions have spatially localized 
structures that were previously believed to only exist on domains with large aspect ra¬ 
tios. We show that one of the newly found steady state solutions underpins the temporal 
intermittencies, i.e., high energy dissipation episodes of the flow. More precisely, it is 
shown that each intermittent episode of a generic turbulent trajectory corresponds to its 
close passage to this equilibrium solution. 

1. Introduction 

The idea that turbulent fluid flow can be studied in terms of the invariant solutions 
of Navier-Stokes equations dates back to 1940s (Hopf 1948). Examples of such invariant 
solutions (also called ‘exact coherent states’ or ‘recurrent flows’ in fluid dynamics liter¬ 
ature) are steady states (or equilibria), traveling waves (or relative equilibria), periodic 
orbits, relative periodic orbits and partially hyperbolic tori. A turbulent flow is then 
viewed as a walk from the neighborhood of one invariant solution to the other. If all 
these solutions are unstable, the turbulent trajectory never settles down and its itinerary 
becomes desperately complex. 

Initially a mathematical endeavor, this view has been put to practice over the last 15 
years in many experimental and numerical studies, providing insights that are beyond 
the reach of purely statistical descriptions of turbulence (see Kawahara et al. (2012) and 
Cvitanovic (2013), for reviews). 

The invariant solutions of Navier-Stokes equations often exhibit complex spatiotempo- 
ral behaviors, and hence analytic expressions for them are unavailable. Their numerical 
computation was first explored in the pioneering works of Marcus & Tuckerman (1987) 
and Nagata (1990). Nagata found the first non-trivial equilibrium and traveling wave 
solutions of the plane Couette flow (Nagata 1990, 1997). Nagata’s approach to finding 
these solutions is the following. Through a hnite truncation, the invariant solutions are 
formulated as the zeros of a large system of nonlinear algebraic equations F(a) = 0. The 
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solution a is then found through Newton-Raphson iterations. Depending on the type 
of the discretization, the finite vector a contains either the Galerkin coefficients or the 
collocation values. 

Nagata’s approach underlies the state-of-the-art methods for finding invariant solutions 
of fluid flow (see, e.g., Waleffe (1998); Kawahara & Kida (2001); Viswanath (2007); 
Gibson et al. (2009); Willis et al. (2013)). Such Newton-type approaches suffer from two 
major drawbacks: 

(i) Obtaining the Newton direction Sa. requires solving the linear system of equations 
VF(a)(5a = —F(a). Even for moderate Reynolds numbers, this linear equation is too 
large to be solved accurately (within numerical precision) in reasonable computation 
time. 

(ii) The convergence of the Newton iterations is not generally guaranteed. The itera¬ 
tions converge only if very good initial guesses are provided. 

As regards item (i), the linear system is typically too large for the matrix VF(a) 
to be even explicitly formed (due to memory considerations). To address this issue, 
the generalized minimal residual (GMRES) method is used to approximate the solution 
of the linear system (Saad & Schultz 1986). The GMRES method replaces the exact 
Newton direction with its least square approximation within a Krylov subspace. Forming 
this subspace only requires the action of the matrix VF(a) on certain vectors, hence 
avoiding the formation of the matrix itself. The resulting method is often referred to as the 
Newton-GMRES iterations. Although computationally feasible, the Newton-GMRES 
iterations are still a formidable numerical undertaking in terms of implementation and 
computation time. 

As mentioned in item (ii) above, even when the exact Newton descent direction 5a is 
known, the iterations are not guaranteed to converge. Newton’s domain of convergence in 
the context of Navier-Stokes equations is typically small, demanding a very good initial 
guess for the iterations to converge (Waleffe 1997). The Newton iterations, therefore, are 
very effective for bifurcation analysis where the invariant solutions bifurcate from known 
solutions (Tuckerman et al. 2000; Waleffe 2003; Kreilos & Eckhardt 2012). In this case, 
the known solution (corresponding to the old bifurcation parameter value) is used as the 
initial guess for the Newton iterations. Even then, the method might fail to determine 
solutions far from the bifurcation point. 

In practice, good initial guesses for Newton iterations are obtained through heuris¬ 
tic methods, such as close recurrences of a turbulent trajectory (Auerbach et al. 1987; 
Kawahara et al. 2006; Viswanath 2007; Cvitanovic & Gibson 2010). Such heuristics lack 
a solid mathematical basis. As a result, many important invariant solutions may and will 
remain undiscovered (see Section 5, below). 

Viswanath (2007) used the hook step together with the Newton-GMRES method which 
significantly improved the convergence of the iterations. The underlying idea of the hook 
step is to approximate the Newton direction within the Krylov subspace with the con¬ 
straint that 5a belongs to a ‘trust region’ such that |5a| is smaller than some prescribed 
value e. The constraint |5a| < e ensures that the linearization VF(a)5a involved in 
Newton iterations is a valid approximation (Dennis Jr. & Schnabel 1996). The result¬ 
ing method, i.e. the Newton-GMRES-hook iterations, is state-of-the-art in computing 
invariant solutions in the context of fluid flow. While improving the convergence chance 
of the Newton-GMRES iterations, the hook step still requires a good initial guess. 

Such drawbacks have slowed down research progress in transitional and moderate 
Reynolds number turbulence. As the study of fluid turbulence through its invariant 
solutions reaches a mature state, it is time to revisit the methods by which these in¬ 
variant solutions are found. Better methods should ideally scale linearly (in the sense 
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of computational complexity) with the number of degrees of freedom and be universally 
convergent. The term universally convergent implies that every initial guess converges to 
some solution. Universally convergent methods, however, do not immediately guarantee 
the identihcation of all possible solutions from a finite set of initial guesses. 

The main goal of the present paper is to propose one such new direction. We develop 
an adjoint-based method for computing steady state and traveling wave solutions of 
incompressible Navier-Stokes equations with periodic boundary conditions. In particular, 
adjoint partial differential equations (PDEs) are derived whose trajectories converge to 
(relative) equilibria of the Navier-Stokes equations. 

Adjoint equations appear naturally in optimal control, where the governing equations 
act as constraints (see Gunzburger (2002), for a survey of applications in fluid dynamics). 
Pringle & Kerswell (2010) and Monokrousos et al. (2011) used adjoint-based optimization 
to find the optimal route to turbulence in transitional shear flows (also see Kerswell et al. 
(2014) for a review). To investigate the regularity of the Naveir-Stokes equations, Ayala 
& Protas (2014) use adjoint-based optimization to find the initial conditions that result 
in maximal palinstrophy growth. 

To the best of our knowledge, adjoint-based methods have not been used to find in¬ 
variant solutions of the Navier-Stokes equations. They have, however, been utilized in 
the context of nonlinear wave equations. For instance, Ambrose & Wilkening (2010&) 
formulate the time-periodic solutions of a unidirectional water wave equation as the 
minima of an appropriate functional. The minima are found iteratively by a ‘steepest 
descent’ method. At each iteration, the descent direction is obtained as the solution of 
a backward-time adjoint PDE (also see Ambrose & Wilkening (2010a)). Yang & Lakoba 
(2007) realize that, for equilibrium solutions, the descent direction can be expressed 
explicitly. Namely, they show that the solitary solutions of nonlinear wave equations 
coincide with the asymptotically stable equilibria of an adjoint PDE. 

The case of Navier-Stokes equations is more involved due to the presence of nonlo¬ 
cal pressure gradients enforcing incompressibility. This calls for a special treatment as 
presented here. We also extend the adjoint-based approach to finding traveling wave 
solutions with a priori unknown wave speeds. 

Recently, Chandler & Kerswell (2013) and Lucas & Kerswell (2014) carried out the 
most exhaustive search for invariant solutions of a chaotic two-dimensional Kolmogorov 
flow using Newton-GMRES-hook iterations. Here, we apply our adjoint-based method to 
the same Kolmogorov flow and discover several new steady state and traveling wave so¬ 
lutions. We show that some of these new solutions underpin the temporal intermittencies 
associated with short-lived, rapid growth of energy dissipation. 

This paper is organized as follows. The preliminary material is reviewed in Section 2, 
in a rather general setting. The explicit form of the adjoint equations for Navier-Stokes 
equations are presented in Section 3. We test the adjoint-based approach in Section 4 
on a two-dimensional Kolmogorov flow and also carry out a thorough comparison with 
Newton-GMRES-hook iterations. In Section 5, the temporal intermittency of the Kol¬ 
mogorov flow is studied in terms of its invariant solutions. Finally, Section 6 contains our 
concluding remarks and an outline of future research directions. 


2. Preliminaries: Newton descent vs. adjoint descent 

In this section, we review Newton and adjoint descent methods in a general framework. 
We find this exposition helpful for understanding the Navier-Stokes adjoint equations 
presented in Section 3. 
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Consider the partial differential equation (PDE) 

dtU = F{u), (2.1) 

where the real valued vector field u(x, t) is a function of space x and time t, and F is 
a nonlinear differential operator acting on an inner product function space R. We seek 
equilibrium (or steady state) solutions of this PDE, i.e., time-independent vector fields 
u = u(x) such that F(u) = 0. 

The finite dimensional counterpart of (2.1) is the system of ordinary differential equa¬ 
tions (ODEs) da/dt = F(a) where a G and F : —)■ While our derivations 

are in the infinite-dimensional setting, we will occasionally invoke this finite-dimensional 
analogue to elaborate the ideas. 

We define the ‘weighted’ inner product 

(u,u ').4 = (u,^u')z, 2 , (2.2) 

for any two functions u, u' G H, where A denotes a real-valued, positive-definite, self- 
adjoint (with respect to the inner product operator. The prime signs should 

not be confused with derivatives; we use them here to distinguish different functions. The 
inner product (2.2) defines a natural norm given by 

l|u|U = V (u,u)^, (2.3) 

for u G H. The choice of A is somewhat arbitrary and should be judiciously chosen for 
a specific problem. For the Navier-Stokes equation, for instance, we will use an operator 
which is closely related to the inverse of the Laplacian (see Section 3). For this exposition, 
however, we set A to identity, which renders (2.2) as the usual inner product. 

With this setting, we seek a second PDE 

9t-u=G(u), (2.4) 

such that its solutions converge to the equilibrium solutions of (2.1) as the fictitious time 
T tends to infinity. To ensure this convergence, the differential operator G is designed in 
a specific way such that 

||F(u(t))||l 2 ^ 0, as T ^ oo, 

where u(t) = u(-, r) denotes a solution of (2.4). Newton and adjoint descents correspond 
to two different choices of the operator G. 

The finite-dimensional analogue of (2.4) is the fictitious-time ODE da/dr = G(a) with 
G : ^ 


2.1. Newton descent 

The main observation here is that discrete Newton iterations are an explicit Euler ap¬ 
proximation of a continuous fictitious-time dynamical system (Saupe 1988). This well- 
known fact is rarely mentioned in the literature, prompting Smale (1981) to write “T/iis 
geometry is based on an idea which is known, but not usually explicated. Namely, New¬ 
ton’s method for solving f{z) = 0 is an Euler approximation to the ordinary differential 
eguation... Here, we first derive the continuous hctitious-time Newton method in the 
infinite-dimensional setting and then show how the discrete Newton iterations follows 
from there. 

The norm ||F(u)||j ;,2 evolves along the trajectories of (2.4) according to 
a.|lF(u)|!i. =2(/:(u;G(u)),F(u))i2, 


(2.5) 



Invariant solutions of the Navier-Stokes equation 


5 


where the linear operator £(u; •) is the Gateaux derivative of F evaluated at the state u, 
and is defined by 


£(u;u') := lim 

£- 5-0 


F(u + eu') -F(u) 


Vu' e H. 


( 2 . 6 ) 


Requiring G to satisfy 

£(u;G(u)) = -F(u), (2.7) 

one obtains the continuous fictitious-time Newton descent in infinite dimensions. Substi¬ 
tuting this expression in (2.5), we obtain 

a.||F(u(r))|ii. =-2|lF(u(r))||i., (2.8) 


which has the exact solution ||F(u(t))|| = l|F(u(0))|| j^ 2 exp(—r). This shows that, 

as long as (2.7) has a solution, ||F(u(T))||i 2 decays to zero exponentially fast along 
the solutions of (2.4). In other words, the continuous fictitious-time Newton method 
converges to an equilibrium solution of (2.1) for almost all initial conditions (Saupe 
1988; Cvitanovic & Lan 2003; Lan & Cvitanovic 2004). 

For a given state u, Eq. (2.7) is a PDE that can in principle be solved for G. In practice, 
it is approximated through some finite truncation to take the form VF(a)G(a) = —F(a), 
with a being the Galerkin coefficients. The solution G of this large, but finite-dimensional, 
linear system is often approximated by some variant of generalized minimal residual 
(GMRES) iterations (Trefethen & Bau 1997). 

Furthermore, the PDE (2.4) is discretized in time to yield the explicit Euler step 


Uj+i = Ui -l-^rG(ui), 


(2.9) 


with 0 < (5t < 1. The standard Newton iterations correspond to St = 1, while damped 
Newton iterations adjust St in order to ensure that the if norm decreases (Boyd & 
Vandenberghe 2004). 

The above discrete iterations do not guarantee the global convergence that the continu¬ 
ous fictitious-time Newton descent does (Saupe 1988). As mentioned in the Introduction, 
the convergence of Newton iterations is only guaranteed if the initial guess is sufficiently 
close to an equilibrium. In fact, if such a close initial guess is available, the convergence 
to the equilibrium is super-exponential (Deuflhard 2011). 

Anecdotal evidence suggests that, in the context of fluid flow, the basin of attraction 
of Newton iterations is rather small (Waleffe 1997; Viswanath 2007). This domain can be 
enlarged by choosing a higher order scheme for temporal discretization of (2.4). This is, 
however, computationally too demanding for large dimensional systems, such as turbulent 
fluid flow. 


2.2. Adjoint descent 

The adjoint descent corresponds to an alternative choice of G that may be expressed 
analytically, thus avoiding the approximation involved in the Newton-GMRES iterations. 
Moreover, the adjoint direction can be evaluated at a relatively low computational cost, 
rendering a higher order time discretization of (2.4) feasible. 

To express the adjoint direction, note that (2.5) can be written as 

5.||F(u)|li. = 2(G(u),£t(u;F(u)))^2, (2.10) 

where the adjoint jC^(u; •) is the linear operator satisfying 

(£(u;u'),u")i2 = (u',£^(u;u"))i2, Vu,u',u" e L^. 


( 2 . 11 ) 
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Requiring 

G(u) = -£t(u;F(u)), (2.12) 

ensures the descent of the norm ||F(u(r))||i 2 along the trajectories of (2.4), since 

9.|lF(u)|ii. =-2||/:t(u;F(u))||2, <0. (2.13) 

We refer to G given by (2.12) as the adjoint direction. The finite-dimensional analogue 
of this adjoint direction is G(a) = — [VF(a)]^ F(a), where T denotes the usual matrix 
transposition. 

The downside of using the adjoint descent is that, unlike the continuous fictitious¬ 
time Newton descent, equation (2.13) does not guarantee an exponential decay of the 
norm. In fact, we observe a relatively slow decay in Section 4.3, and propose a hybrid 
adjoint-Newton algorithm to rectify this drawback. 

In Section 3, we derive the explicit form of the adjoint operator for the incompressible 
Navier-Stokes equations. 


3. Adjoint descent for the Navier—Stokes eqnations 

Consider the incompressible Navier-Stokes equation in non-dimensional variables, 

dju = —u • Vu — Vp-I--I-f, (3.1o) 

V-u = 0, (3.16) 

defined on the spatial domain x = {xi,X2, x^) e I? = [ 0 , Li] x [ 0 , L 2 ] x [ 0 , L 3 ] with periodic 
boundary conditions. The parameter v is the inverse of the Reynolds number, v = Re~^. 
We assume that the time-independent forcing term f = f (x) is divergence-free, V • f = 0. 

The goal is to find an adjoint PDF such that, along its solutions (u(r),p(T)), the 
norm of the right-hand-side of (3.1a) decays to zero monotonically, while ensuring that 
V - u = 0 for all r. We derive the adjoint equations for the general norms of the type (2.3). 
In Section (3.3), we make a specific choice for the norm that is used in the subsequent 
computations. 


3.1. The adjoint descent equation for equilibria 

Define 

Fo(u) = — u • Vu — Vp -b i/Au -I- f, (3.2) 

which is the right-hand-side of Eq. (3.1a). The reason for using the subscript 0 in Fq 
becomes clear in Section 3.2, where we derive the adjoint equation for traveling waves. 
We seek steady states u such that Fo(u) = 0. 

To this end, we derive an adjoint PDE, such that, along its trajectories u(t), the norm 
||Fo(u(t))||_4 decreases monotonically and converges to zero asymptotically. Our choice of 
the norm does not compromise the accuracy of the resulting equilibrium solutions of the 
Navier-Stokes equations: due to the positive-definiteness of the operator A, |jFo(u)m = 
0 if and only if Fq (u) = 0 

Leaving the derivation to Appendix A, the adjoint descent equations read 

= - { [Vu" -b (Vu")^] u - Vp" + i/Au"} , (3.3a) 

V • u" = 0, V • u = 0, (3.36) 

where u" = Fo(u), i.e., 

u" = —u • Vu — Vp -b :^Au + f, (3.4) 
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and we use the shorthand notation u" := ^u". 

Two pressure-type terms p and p" appear in the adjoint equation. As in the case of 
the Navier-Stokes equation, the pressure terms play the role of the Lagrange multipliers, 
enforcing the divergence-free constraints on u and u" (cf. Appendix A). Taking the 
divergence of equation (3.4) we get 

Ap = —V • [u • Vu], (3.5) 

enforcing V • u" = 0. Similarly, taking the divergence of (3.3a) yields 

Ap" = V • { [Va" + (Vfi")T] u} , (3.6) 

ensuring V • u = 0. 

Due to the periodic boundary conditions, these divergence-free constraints can be easily 
enforced by a projection operator in the spectral space. In Appendix C, we derive the 
Fourier space representation of the adjoint equations (3.3). This spectral representation 
closely resembles that of the Navier-Stokes equations. Therefore, an existing pseudo- 
spectral code for the Navier-Stokes equations can also be used, with straightforward 
modihcations, to solve the adjoint equations. Furthermore, the adjoint equations (3.3) 
enjoy the same degree of parallelizability as the Navier-Stokes equations. 


3.2. Adjoint descent for the traveling waves 

In the absence of forcing, the Navier-Stokes equation (3.1) is invariant under Galilean 
translations. This symmetry allows for the existence of traveling wave (or relative equi¬ 
librium) solutions. A traveling wave has the general form u(x,t) = u(x — ct), where 
c = (ci,C 2 ,C 3 ) S is a constant wave velocity. The forcing term f may break the 
translational symmetry in one or more directions, in which case, the wave speed Ci cor¬ 
responding to the symmetry-broken direction Xi is identically zero. 

Substituting the traveling wave solution u(x—ct) into the Navier-Stokes equation (3.1), 
gives the time-independent equations 

— u • Vu — Vp-I--I-f-f c • Vu = 0, V • u = 0. (3.7) 

Therefore, finding traveling wave solutions to the Navier-Stokes equation amounts to 
finding the kernel of the operator 

Fc(u) = —u • Vu — Vp -I- vAu. -|- f -I- c • Vu. (3.8) 

If the wave velocity c is known, the solutions to Fc(u) = 0 can be found through 
an adjoint descent equation similar to (3.3). In general, however, the wave velocity c is 
unknown. To address this general case, we allow c to be a function of the fictitious time 
r, and enforce its derivative c = dc/dr to change in such a way that ||Fc(T-)(u(T))m 
decreases monotonically to zero as t increases. 

Leaving the details to Appendix B, we obtain the set of adjoint descent equations 

a^u = - { [Vfi" (Vii")T] u - Vp" -f vAxjf'] + c • Vii", (3.9a) 


V • u" = 0, V • u = 0, 

where u" = Fc(u), i.e.. 


(3.96) 

(3.9c) 


u" = —u • Vu — Vp -I- i^Au -I- f -I- c • Vu, (3.10) 

and u" = Au". The pressure-type functions p and p" satisfy equations (3.5) and (3.6), 
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respectively, with u" given by (3.10). The spectral representation of the adjoint equa¬ 
tions (3.9) is given in Appendix C. 

The adjoint descent for traveling waves includes the evolution equation (3.9&) for the 
unknown wave velocity c. If the symmetry in Xi direction is broken due to the forcing 
term f, the corresponding wave speed Ci is identically zero and the equation for Ci is 
eliminated from (3.96). 


3.3. Choice of the norm 

The conventional choice for the norm || • ||_4 is the norm, corresponding to A being the 
identity operation. We find, however, that the resulting equations are stiff, requiring very 
small time-steps for their numerical integration. This is in line with a similar observation 
made by Yang & Lakoba (2007) in the context of nonlinear wave equations. Here we 
choose the H~^ norm that renders the adjoint PDEs less stiff by damping the higher 
Fourier modes. 

This choice is somewhat arbitrary, but we find it to work very well for the Kolmogorov 
flow presented in Section 4.1. Moreover, the implementation of the resulting equations 
in the Fourier space is straightforward. 

The H~^ can be computed as follows. Let q = (u,p) be the velocity-pressure pair and 
define the operator A through its action in the Fourier space, 

( 3 . 11 , 


where q denotes the Fourier transform of q and k = (27rfci/Li, 2Trk2/L2, 2TTk2,/L^), with 
ki € being the wavenumber. Note that the operator A commutes with the divergence 
operation, i.e., A(V • u) = V • Au = 0. 

With this choice and using Parseval’s identity, the inner product (2.2) can be written 
explicitly in the Fourier space as 




= E 


q(k) • q'(k) 

l + |k|2 


(3.12) 


where * denotes complex conjugation. This inner product induces the H ^ norm ||q||^_i = 

(q, q)ff-i- 

The operator A appears in the adjoint equations (3.3) and (3.9) through u" = Au", 
which can be readily computed in the Fourier space through Fq. (3.11). 


4. Equilibria and traveling waves for Kolmogorov flow 

We test the performance of the adjoint descent equations on a two-dimensional Kol¬ 
mogorov flow. This flow (with identical geometry and parameters as the ones chosen here) 
was studied thoroughly by Chandler & Kerswell (2013), providing a basis for comparison 
with our results. 

In Section 4.3, we introduce the hybrid adjoint-Newton method which requires the 
Newton-GMRES-hook (or NGh, for short) iterations. The NGh iterations are explained 
in detail by Ghandler & Kerswell (2013). Our implementation and choice of parameters 
are identical to theirs with the exception that we use the velocity formulations of the 
Kolmogorov flow, while Ghandler & Kerswell (2013) use the vorticity formulation. We 
chose to use the velocity formulation since the adjoint equations (3.3) and (3.9) are 
derived for the general three-dimensional flow, and hence are in the velocity-pressure 
form. 




9 


Invariant solutions of the Navier-Stokes equation 
4.1. Kolmogorov flow 

Considered the forced Navier-Stokes equation (3.1) on the two-dimensional torus x = 
{xi,X 2 ) S = [0, 27r] X [0, 27r] with the forcing f(x) = sin(na: 2 )ei, where ei = (1,0)^ 
and n is a positive integer. This choice of forcing yields the two-dimensional Kolmogorov 
equation 

dtu = —u • Vu — Vp -1- -I- sin(na; 2 )ei, (4.1a) 


V-u = 0, (4.1&) 

where, u = {ui,U 2 ) and u = 1/Re with Re being the Reynolds number. 

Here, we review some of the relevant properties of this equation that can also be found 
in Platt et al. (1991) and Chandler & Kerswell (2013). The Kolmogorov equation has an 
equilibrium solution that can be expressed analytically as 

Re 

ui{xi,X2) = ^sm.{nx2), U2{xi,X2) = (4.2) 

at any Reynolds number Re. Following Chandler & Kerswell (2013), we refer to this 
solution as the laminar state Eq. 

The instantaneous energy E, energy dissipation D and energy input / of a state u are 
defined, respectively, by 


m 

Dlt) 

I(t) 


2(2 




1 


1 


[[ Ui{xi,X2,t)sm{nx2)dx, 
(2TTy jJy2 


JJ 2 


(4.3a) 

(4.36) 

(4.3c) 


where | • | denotes the usual Euclidean norm and uj is the non-zero component of the 
vorticity field V x u. One can readily show that 


dE 

dt 


= I-D, 


(4.4) 


along an arbitrary solution of the Navier-Stokes equation. For equilibrium and traveling 
wave solutions, the energy is time-independent and therefore I = D = const. For the 
laminar solution Eq, for instance, we have 


Elam — 


Re^ 


Re 


(4.5) 


For n = 1, the laminar state is the global attractor for all Reynolds numbers, and 
therefore all solutions converge to it asymptotically (see, e.g., Marchioro (1986) and 
Foias et al. (2001) (Section 3.1)). For n > 1 and large enough Reynolds numbers, the 
laminar state becomes unstable. The numerical study of Platt et al. (1991) suggests that 
for n = 4 and large enough Reynolds numbers, all invariant solutions of Kolmogorov 
equation become unstable, resulting in a chaotic attractor. More recently, Chandler & 
Kerswell (2013) confirmed this observation. To be in agreement with these studies, we 
also choose the forcing wavenumber n = 4. Our numerical results are carried out at 
Re = 40, unless otherwise is mentioned explicitly. Chandler & Kerswell (2013) also carry 
out their most exhaustive search for invariant solutions at Re = 40. 
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Figure 1. Typical energy spectrum of a turbulent trajectory at Re = 40. 


4.2. Temporal and spatial discretization 

To evaluate the right-hand-side of the Navier-Stokes equation (4.1) and the adjoint 
equations (3.3) and (3.9), we use a standard pseudo-spectral method with 2/3 dealiasing. 
Chandler & Kerswell (2013) consider Re = 40, 60, 80 and 100 using a uniform 256 x 256 
Fourier modes for all Reynolds numbers. For Re = 40, this resolution is unnecessarily 
high. Instead, we use 128 x 128 Fourier modes for Re = 40. The resulting energy spectra 
(cf. Fig. 1), and the fact that we have been able to reproduce the (relative) equilibria 
found by Chandler & Kerswell (2013), reassures that 128^ modes are sufficient. 

For the time integration of the adjoint equations (3.3) and (3.9), one can take an explicit 
adaptive Euler time-step. Denoting the right-hand-side of either adjoint equation by G, 
the Tth Euler step reads 

Ui+i = Uj -I- Sn G(ui), 

where the length Sri of the time step is chosen small enough such that the residue 
||F(uj+i)|U is less than ||F(uj)|U. 

Yang & Lakoba (2007) derive an upper bound for the admissible time step Com¬ 
puting this upper bound for Navier-Stokes equation is, however, not straightforward. In 
practice, therefore, one starts with a large value Sri and decreases it incrementally until 
the residue ||F(ui+i )||_4 reduces compared to the previous iteration. 

To gain accuracy, however, we use a higher order numerical scheme for our time in¬ 
tegrations. This is computationally feasible due to the low cost of evaluating the right- 
had-sides of equations (3.3) and (3.9) by the pseudo-spectral method. More specifically, 
we use the embedded Runge-Kutta scheme RK5(4) of Dormand & Prince (1980). This 
scheme allows for an automatic adaptive time step-size. Roughly speaking, RK5(4) takes 
forth and hfth order Runge-Kutta (RK) steps. The fifth order is eventually used for 
the time stepping. The forth order prediction is used to adjust the step-size as follows. 
Let err denote the absolute difference between the forth and the fifth order predictions. 
Then the step-size is chosen such that err < atol + |ui|rtol. The prescribed constants 
atol and rtol are the absolute and relative errors, respectively. We refer the reader 
to Press et al. (2007) (Section 17.2) for further details and an implementation of the 
RK5(4) scheme. This integrator is also implemented in the MATLAB function ode45. 
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Figure 2. (a) The vorticity field for the initial condition u(x, 0) = (cos(2a;2), cos( 2 :i))^. (b) The 
initial condition evolved to time r = 500 under the adjoint equation (3.3) to obtain u(x, 500). (c) 
The result u(x, 500) shown in panel (b) is used as the initial guess for the Newton-GMRES-hook 
iterations which converged after 7 iterations with residue 5.98 x 10“^^. This equilibrium solution 
is labeled in Table 1 as E 4 . (d) The evolution of the L^-residue ||Fo(u(-,T))||i 2 as the initial 
condition in (a) evolves under the adjoint equation from r = 0 to r = 500 


For integrating the adjoint equations, we choose atol = rtol = 10“^°. Even with this 
conservative choice, step sizes as large as 10 fictitious-time units were taken by RK5 (4). 

We also use RK5(4) for the temporal discretization of the Navier-Stokes equation (4.1). 
However, as the time step sizes were much smaller for integrating this equation, we used 
the less stringent error tolerances atol = rtol = 10“®. 

4.3. Hybrid adjoint-Newton iterations 

We find that the adjoint equation (3.3) does take arbitrary initial conditions to Navier- 
Stokes equilibria. The convergence is, however, rather slow. To demonstrate this, we take 
the initial condition u(x, 0) = (cos(2a:2), cos(a:i)) (see Fig. 2(a)) and evolve it under the 
adjoint equation (3.3) to time r = 500. The result is shown in Fig. 2(b). This integration 
took 54 seconds (on an iMac with a single processor: Intel Core i5, 2.9 GHz). The Lf 
norm of the residue Fq, however, decreases to ~ 5 x 10“^ approximately, indicating that 
much longer integration times are required to reduce the residue sufficiently enough, say 
to 10-1°. 

This slow convergence is due to nearly marginal, a priori unknown eigenmodes of the 
adjoint operator (Lakoba & Yang 2007). For particular wave equations, ad hoc methods 
have been proposed to eliminate these modes, and hence speed up the convergence (Yang 
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Figure 3. State space depiction of the hybrid adjoint-Newton method. An initial guess may be 
too far away from the desired equilibrium solution EQ for the Newton-GMRES-hook iterations 
(NGh) to converge to it. The adjoint descent may eventually converge to EQ, but alone it would 
take a long integration time to do so. In the hybrid adjoint-Newton method the adjoint descent 
takes the initial guess to the domain of attraction of NGh. Once there, NGh converges to EQ 
in a few iterations. 


& Lakoba 2007). Due to the complexity of the Navier-Stokes equations, it is unclear how 
such mode elimination techniques could be employed here. 

Using the state u(x, 0) = (cos(2a:2), cos(xi)) directly as the initial guess for the NGh 
iterations does not converge to an equilibrium either: after the first 20 iterations, the 
residue plateaued around 2 x 10“^ and remained so for the 75 iterations that were carried 
out. 

Instead, when we use the result of the adjoint integration, i.e. u(x, r = 500), as the 
initial guess for NGh iterations, it converges after seven iterations with residue 5.98 x 
10“^^ (Fig. 2(c)). These seven iterations took 156.84 seconds. 

This suggests that, although the convergence of the adjoint equation to the equilibrium 
is slow, it takes generic initial guesses to the vicinity of the equilibrium at a relatively 
low computational cost. This can be seen by comparing panels (b) and (c) of Fig. (2), 
which shows that the state u(x, r = 500) is quite similar to the equilibrium and thus 
likely to lie in its NGh domain of attraction. Indeed, switching to NGh took this state 
to the equilibrium solution within a few iterations. 

Based on this observation, we propose the following hybrid adjoint-Newton iterations, 
sketched in Fig. 3. We take a prescribed positive real number tq and a positive integer N. 
An initial condition u(x, 0) is then integrated under the adjoint equation (3.3) to obtain 
u(x. To). This state is then fed into the NGh algorithm and N iterations of NGh are 
performed. The output is then fed back into the adjoint equation as the initial condition. 
This hybrid loop is repeated until the residue |jFo||L 2 decreases below some prescribed 
tolerance tol. This procedure is summarized in Algorithm I. The hybrid adjoint-Newton 
iterations for finding traveling waves are similar except that the adjoint equations (3.9) 
are solved at each iteration (see Algorithm 2). 

For the computations reported below, we set tq = 100, tol = 10“^° and N = 1. Since 
NGh steps are relatively expensive, we only take one NGh step (N = 1) per iteration 
of adjoint-Newton. At Reynolds number Re = 40, for instance, the integration of the 
adjoint equations to r = 100 took approximately 10 seconds on average while each NGh 
step took approximately 50 seconds on average. 
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Algorithm 1: Hybrid adjoint-Newton algorithm for finding equilibrium solutions of 
the Navier-Stokes equation. 

Input: Tojtol G IR+, NG N, state Ug 
while ||Fo(uo)||l 2 > tol do 

Integrate the adjoint equation (3.3) for tq fictitious-time units with the initial 
condition u(-,0) = Ug 

Ug <- u(-,Tg) 

if ||Fo(ug)||i2 < tol then 

L STOP 

for (/c = 1 to N) do 

Take one Newton-GMRES-hook step to get NGh(ug) 

Ug ^ NGh(Ug) 
if ||Fo(ug)|ji2 < tol then 

[ L STOP 

Output: Ug 


Algorithm 2: Hybrid adjoint-Newton algorithm for finding traveling wave solutions 
of the Navier-Stokes equation. 

Input: Tg,tol G K+, NG N, state ug, wave speed cg 
while ||Fco(ug)||L 2 > tol do 

Integrate the adjoint equation (3.9) for rg fictitious-time units with the initial 
conditions u(-, 0) = Ug and c(0) = Cg 

Ug ^ u(-,Tg) 

Cg ^ c(Tg) 

if l|Fco(ug)||L 2 < tol then 
L STOP 

for (A: = 1 to N) do 

Take one Newton-GMRES-hook step to get NGh(ug,Cg) 

(ug,Cg) <- NGh(ug,Cg) 

if ||Fco(ug)||i2 < tol then 

L L STOP 

Output: Ug,Cg 


4.4. Equilibrium solutions 

In this section, we report the equilibria found by the hybrid adjoint-Newton iterations 
(Algorithm 1). We also compare its performance with that of pure NGh iterations without 
the adjoint step. 

To this end, we consider the set of initial guesses 

Ug’”!’™^) = (cos(TO2a;2),cos(mia:i)), (4.6) 

for a range of integers mi and m 2 . We find that for m 2 > 4, both Newton-GMRES- 
hook and the hybrid adjoint-Newton iterations converge to the laminar equilibrium Eq. 
Therefore, we restrict the range of the integers to 1 < mi, m 2 < 4. This leads to 16 
distinct initial guesses. 

These initial conditions are divergence-free by construction, consistent with our incom¬ 
pressible Navier-Stokes setting. They are also explicit, rendering the following results re¬ 
producible. Furthermore, the generic nature of the initial conditions illustrates the main 
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Figure 4. The decay of the residue (i.e., the L^-norm of the right hand side) using 
Newton-GMRES-hook iterations (a) and the hybrid adjoint-Newton iterations (b). 



(3,4) 


Figure 5. The convergence diagram. Double indices (mi,m 2 ) label the initial guesses (4.6). 
A dashed black arrow indicates that both the hybrid adjoint-Newton and the NGh iterations 
converge to the same equilibrium. A dotted blue arrow indicates that the NGh iterations converge 
to a different equilibrium from the one reached by the hybrid adjoint-Newton (red arrow). A 
solid red arrow (with no blue arrow originating from the same initial guess) indicates that only 
the hybrid adjoint-Newton converged to the equilibrium, i.e., the NGh iterations alone failed to 
converge. 


advantage of the adjoint descent method, namely that its convergence does not require 
good initial guesses. For brevity, we shall refer to the states by their indices 

(mi,TO2). 

Using the hybrid adjoint-Newton iterations for equilibria (see Algorithm 1), all 16 
runs converged, resulting in 10 distinct equilibria. These equilibria are listed as Eq to 
Eg in Table 1, where Eq is the laminar state (4.2). All initial guesses converged to an 
equilibrium within the first 10 iterations of hybrid adjoint-Newton iterations, except 
initial guess (2,1) that took 17 iterations (see Fig. 4). Fig. 5 shows the convergence 
diagram, connecting each initial guess (mi, mg) to the resulting equilibrium solution. 

Using NGh alone, only 9 out of 16 runs converged. They converged to 6 distinct equi¬ 
libria: Eg, El, Eg, Es, Eg, Eig and Eu (see Fig. 5). This comes initially as a surprise 
since the exhaustive search carried out by Chandler & Kerswell (2013) only returned 
a single equilibrium (i.e., Ei). This can be explained, however, in terms of the method 
used for initiating the NGh iterations. Chandler & Kerswell (2013) use recurrences to 
find initial guesses for their searches, as opposed to the generic initial guesses used here. 
Recurrences only happen in a subset of the state space where a generic turbulent tra¬ 
jectory spends most of its lifetime. As a result, using recurrences for initiating the NGh 
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Equilibrium I=D 

E 


uJi dim E’ 

Eo 

1.25000 

1.56250 

2.35340 

0.0 

38 

El 

0.12732 

0.76108 

0.24929 2.56010 

9 

E 2 

0.15406 

0.47313 

0.66722 

0.0 

13 

E 3 

0.26287 

0.45173 

0.84744 

0.0 

27 

Ei 

0.08433 

0.57317 

0.62697 

0.0 

5 

Es 

0.26661 

0.55183 

0.30119 

0.0 

19 

Ee 

0.26227 

0.45574 

0.86924 

0.0 

22 

E7 

0.07530 

0.58104 

0.58318 

0.0 

5 

Eg 

0.17452 

0.53493 

0.61189 0.00675 

17 

Eg 

0.15315 

0.47396 

0.67006 

0.0 

17 

Eio 

0.61437 

0.91740 

0.85655 : 

1.75030 

30 

Ell 

0.48020 

0.71284 

0.81175 

0.0 

29 

El 2 

0.31049 

0.89032 

0.68928 

0.0 

13 

El 3 

0.26151 

0.50921 

0.42190 : 

1.49660 

16 

El 4 

0.31152 

0.54066 

0.71722 

0.0 

21 

El 5 

0.27070 

0.76570 

0.87760 

0.0 

14 

El 6 

0.34954 

0.61168 

0.83092 0.02382 

24 


Table 1. List of equilibrium solutions at Re = 40. Energy E, energy dissipation D and energy 
input I are defined in (4.3). The leading unstable eigenvalue of the equilibrium is /ii + iwi. The 
dimension of the linear unstable manifold of the equilibrium is denoted by dimiJ“. 



Figure 6 . Vorticity fields of the equilibrium solutions E 5 (left), E \2 (middle) and E 15 (right). 
All panels show the entire domain [0, 27r] x [0, 27r] 


searches might preferentially yield the equilibria ‘close’ to this region. We will return to 
this subject in Section 5, where the temporal intermittency is discussed. 

We also searched for further equilibria using states of the form (sin(TO2a;2), cos{miXi)) 
as the initial guesses for the hybrid adjoint-Newton iterations. This resulted into five more 
equilibria: E 12 to Eiq in Table 1. While all of these additional searches did converge, most 
of them re-converged to previously discovered equilibria, including Eiq and En that were 
only found by NGh iterations when initial guesses (4.6) were used. Figure 6 shows the 
vorticity field of three select equilibrium solutions. 

Some of the equilibria (e.g., £^15 shown in Fig. 6 ) exhibit vertical bands of localized 
vorticity that are separated by an almost zero vorticity background. Such spatially local¬ 
ized equilibria of the Kolmogorov flow were only observed previously on domains with 




























Figure 7. Vorticity fields of the equilibrium solutions at Reynolds number Re = 60 (left), 
Re = 80 (middle) and Re = 100 (right). All panels show the entire [0, 27r] x [0, 27r] domain. 


small aspect ratio a = L 2 ILi (Lucas & Kerswell 2014). The fact that they also exist on 
a domain with aspect ratio a = 1 comes as a surprise. 

Although our focus here is on Reynolds number Re = 40, Fig. 7 showcases select 
equilibria at Re = 60, 80 and 100. These equilibria are computed using the higher 
resolution of 256 x 256 Fourier modes. They are found by our hybrid adjoint-Newton 
method while the previous studies using Newton-GMRES-hook iterations had not been 
able to discover them (Chandler & Kerswell 2013). 


4.5. Traveling wave solutions 

The forcing term in the Kolmogorov equation (4.1) breaks the continuous symmetry in 
the a; 2 -direction. Therefore, only traveling wave solutions of the type u(x — ct) with 
c = (c, 0) are permitted. This reduces the wave velocity equations (3.9&) to the scalar 
equation 


dc 

dr 



(4.7) 


Similarly, the term c • Vu" in (3.9a) reduces to cdxiU” and the term c • Vu in (3.10) 
reduces to cOx^u. 

We search for traveling waves using Algorithm 2 and generic initial guesses discussed 
in the previous section. For the initial wave speed, we used c(0) = 1. Other values of c(0) 
yielded similar results. 

Some of our searches for traveling waves, converged to equilibrium solutions. This is to 
be expected as equilibria are degenerate traveling waves with wave speed c = 0. In fact, 
the adjoint equation (3.9) admits such solutions. This can be readily verified by letting 
u(x) to be an equilibrium solution of the Navier-Stokes equation and setting c = 0. Then 
u" in (3.10) is zero, resulting in vanishing right-hand-sides in equations (3.9). 

Nonetheless, our hybrid adjoint-Newton searches led to 9 distinct traveling waves listed 
in Table 2. Only traveling waves Ti and T 2 had been discovered previously (Chandler & 
Kerswell 2013). Figure 8 shows the vorticity field for three select traveling waves. As in 
the case of equilibria, we find that some of the traveling waves (e.g., Tg in Fig. 8) exhibit 
localized spatial structures, although the domain aspect ration is one. 
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Traveling wave 

c 

I=D 

E 

hi 

LOl 

dim A’ 

Ti 

0.01978 

0.08873 

0.69747 

0.06815 

0.35451 

4 

Tz 

0.00944 

0.08680 

0.63969 

0.45288 

0.02134 

4 

T3 

0.01826 

0.13432 

0.38056 

0.49301 

0.0 

10 

T4 

0.03062 

0.31453 

0.52793 

0.55183 

0.0 

21 

Ts 

0.05266 

0.40027 

0.61651 

0.82112 

0.0 

17 

n 

0.04223 

0.32063 

0.52963 

0.54292 

0.0 

18 

Tt 

0.00208 

0.08058 

0.61200 

0.50598 

0.0 

2 

n 

0.00156 

0.09641 

0.59883 

0.63357 

0.0 

7 

Tq 

0.00642 

0.08867 

0.61764 

0.57681 

0.0 

4 


Table 2. List of traveling wave solutions at Re = 40. The constant c denotes the wave speed. 
Energy E, energy dissipation D and energy input I are defined in (4.3). The leading unstable 
eigenvalue of the traveling wave is + ioJi. The dimension of the linear unstable manifold of 
the traveling wave is denoted by dimi5“. 



Figure 8. Vorticity field for the traveling wave solutions T 4 (left), T7 (middle) and Tg (right). 
All panels show the entire domain [0, 27r] x [0,2^]. 


5. Temporal intermittency in Kolmogorov flow 

In this section, we explore the significance of the invariant solutions, found in the 
previous section, on the global dynamics of the Kolmogorov flow. 

Figure 9 shows the energy input I versus the energy dissipation D for a generic tur¬ 
bulent trajectory computed for 10^ time units and recorded every 0.1 time units. The 
energy input and dissipation of each equilibrium and traveling wave are marked by circles 
and squares, respectively. As mentioned earlier, the energy input and dissipation coincide 
for these invariant solutions, locating them on the diagonal I = D. 

The equilibria and traveling waves assume a wide range of energy input and dissipation. 
The turbulent trajectory, on the other hand, mostly resides in a relatively low energy 
input/dissipation regime. In particular, approximately 85% of this trajectory belongs to 
the I/ham < 0.12 and D/Diam < 0.12 regime, marked by the green square in Fig. 9. For 
the lack of a better term, we refer to this regime as the ‘ergodic sea’. 

At the same time, the turbulent trajectory also experiences sporadic episodes of high 
energy input and dissipation. Such rare, extreme events are usually referred to as tem¬ 
poral intermittency and are ubiquitous in turbulent fluid flow (see, e.g., Batchelor & 
Townsend (1949); Sreenivasan & Antonia (1997)). Non-Gaussian probability distribu- 
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tion of turbulent quantities are a footprint of intermittency that produces the heavy tails 
of the distribution functions (Frisch 1996; Mininni & Pouquet 2010). 

The short lifetime of the intermittent bursts is better seen in Fig. 10 (left panel) where 
the normalized energy input //Iiam is shown as a function of time. The time series for the 
energy dissipation (not shown here) is very similar, except that the intermittent bursts 
of the dissipation occur with a short delay of 0.5 to 1.5 time units relative to the energy 
input bursts. This suggests that, once in a while, the turbulent velocity field mostly aligns 
with the external forcing sin(na; 2 )ei resulting in the growth of the energy input which 
kicks the trajectory out of its ergodic sea. After a short time delay, the energy dissipation 
also increases, bringing the trajectory back to the ergodic regime. 

Spatial intermittency is also a characteristic of turbulent fluid flow which refers to un¬ 
usually large velocity (or vorticity) amplitudes occurring in a relatively small subset of the 
physical domain (see, e.g., Kuo & Corrsin (1971); Schneider et al. (2004); Ruppert-Felsot 
et al. (2009)). Although temporal and spatial intermittencies are sometimes conflated in 
the literature, the relation between the two is not well-understood (Gibbon & Doering 
2003). For our turbulent trajectory, in fact, an appreciable correlation between them was 
not found. Figure 10(b), for instance, shows that the normalized maximum vorticity am¬ 
plitude oscillates rapidly throughout the simulation time, exhibiting no clear correlation 
with the energy input. 

Focusing on the temporal intermittency, we first review the dynamical systems per¬ 
spective of this phenomena. 

5.1. A dynamical systems perspective on temporal intermittency 
Although the Navier-Stokes equations generate an infinite-dimensional dynamical sys¬ 
tem, it is believed that due to the dissipative term z/Au, its solutions converge expo¬ 
nentially fast to a finite-dimensional, invariant set, usually referred to as the inertial 
manifold (see Constantin et al. (1989) for the rigorous definition). The existence of the 
inertial manifold for the Navier-Stokes equation, in its most general form, is an open 
mathematical problem. In practice, however, its existence is often assumed. In fact, this 
assumption underlies the finite Galerkin truncations used in computations (Foias et al. 
1988). 

Some (relative) equilibria and (relative) periodic orbits and portions of their stable and 
unstable manifolds belong to the inertial manifold. For large enough Reynolds numbers, 
these invariant solutions are typically unstable. A generic turbulent trajectory visits the 
neighborhood of an invariant solution for a finite time before it is repelled along its 
unstable manifold towards the neighborhood of another invariant solution (Ruelle 1991; 
Halcrow et al. 2009). This process continues indefinitely in a somewhat unpredictable 
fashion, thereby causing the complex temporal behavior of turbulent trajectories (see 
Fig. 11 for an illustration). 

From this perspective, intermittent episodes are viewed as close passages of the tur¬ 
bulent trajectory to invariant solutions that reside in a ‘less accessible’ region of the 
inertial manifold or, more precisely, the attractor (Holmes 1993). As depicted in Fig. 11, 
such passages are viable along the heteroclinic connections between the invariant solu¬ 
tions (Holmes & Stone 1992; Holmes et al. 2012). 

To characterize the less accessible regions of the attractor, one naturally needs to an¬ 
swer the following question: How frequently is an invariant solution visited by generic 
turbulent trajectories? There is no straightforward, a priori answer to this question (Gvi- 
tanovic et al. 2014, Ghapter 23). There are, however, some characteristics of the invariant 
solutions that are relevant. In Tables I and 2, for instance, we report the dimension of the 
linear unstable manifold (i.e., dimi?“) of each invariant solution. The invariant solutions 
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I / Ilam 

Figure 9. Re = 40. Gray: Turbulent trajectory spanning 10^ time units. Red circles: equilibria. 
Blue squares: traveling waves. The green square marks the region where I/Iiam < 0.12 and 
D/Diam < 0.12. The turbulent trajectory spends 86.62% out of the total 10^ time units inside 
this region. The diagonal / = D is marked by the dashed black line. Equilibria and traveling 
waves with 1/Iiam = D/Diam > 0.32 are not shown. 


El, E 4 , Ej, Ti, T 2 , T 3 , T 7 , Tg and Tg that reside close to the ergodic sea have at most 10 
linearly unstable eigenmodes. The remaining invariant solutions have at least 13 unstable 
eigenmodes and seem to reside further away from the ergodic sea. 

In Tables 1 and 2, we also report the stability exponent fii + iuji of the most unstable 
eigenmode of each invariant solution. It is tempting to assert that the invariant solutions 
with larger fii are less likely to be visited by a generic turbulent trajectory. This is, 
however, not the case. For instance, we have qti — 0.62697 for equilibrium E 4 and fii = 
0.61189 for equilibrium Eg- Equilibrium E 4 is located in the heart of the ergodic sea 
and is, in fact, visited by the turbulent trajectory quite often. Equilibrium Eg, in spite 
of having a similar stability exponent, is rarely visited by the turbulent trajectory. The 
crucial difference between these two equilibria is the dimension of the linear unstable 
manifold which are dimE“ = 5 for E 4 and dimE“ = 17 for Eg. 

In retrospect, the lack of correlation between the stability exponent and the frequency 
of visitations by the turbulent trajectory is to be expected. The stability exponent of 
an invariant solution is a local quantity. As such, for it to be meaningful, the turbulent 
trajectory should already be in the vicinity of the invariant solution. Once there, the 
stability exponent /ii determines how quickly the trajectory will leave the neighborhood. 

At any rate, neither the dimension of the unstable manifold nor the stability exponents 
of an invariant solution decisively determine the frequency at which its neighborhood is 
visited by a turbulent trajectory. Therefore, we take a more direct approach to quantify 
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Figure 10. (a) The energy input 7 as a function of time normalized by the energy in¬ 
put of the laminar state ham = 1.25. (b) The spatial maximum of vorticity magni¬ 
tude, i.e. maXxg'T^ |w(x,t)|, normalized by its value corresponding to the laminar state, i.e., 

max^g.r2 |a;;am(x)| = 10. 



M 


Figure 11 . An initial condition (black square) decays rapidly to the inertial manifold M where 
the dynamics is governed by the (relative) equilibria, (relative) periodic orbits, and portions of 
their stable/unstable manifolds that lie within the inertial manifold. The green dots represent 
equilibria. Highly unstable invariant solutions (e.g. the equilbrium E) are rarely visited by a 
generic trajectory. 


intermittency. Namely, we measure the LZ distance between the turbulent trajectory and 
the computed invariant solutions. 

The Kolmogorov equation is equivariant under a one-parameter family of continuous 
symmetries and a number of discrete symmetries (Sirovich 1987). This implies that each 
solution u(t) has infinitely many equivalent copies. Therefore, when measuring the ‘dis¬ 
tance’ between two states one needs to make an informed choice among the equivalent 
copies of each state. This necessitates a discussion on the symmetries of the Kolmogorov 
flow. 
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5.2. Symmetries of Kolmogorov flow 

The Kolmogorov equation (4.1) is equivariant with respect to 4n discrete symmetries 
and a continuous translational symmetry (Sirovich 1987). We denote the complete set of 
such symmetries by 


{Tiu) {xi,X2) = u(a;i + i, X 2 ) 

( 7 ^u) {xi,X2) = -u(-a;i, -0:2) 


( 5 ™u) {XI,X2) 


({-l)'^ui{{-l)'^xi,X 2 + m-K/n)\ 
U 2 {{-l)'^xx,X 2 + rmT/n)) ' 


(5.1a) 

(5.1fe) 

(5.1c) 


where m £ {0,1,-- - ,2n — 1}. Here, 7^ denotes £-shift in the Xi-direction, TZ denotes 
rotation through tt and S denotes a simultaneous ( 7 r/n)-shift in the a:: 2 -direction and a 
reflection in the xi-direction. 

One can readily verify that the above symmetry operations act on the vorticity field 
according to the following rules: 


(Teuj) {xi,X 2 ) = u}{xi + £,X 2 ) (5.2a) 

(TZuj) {xi,X2 ) = u!{—xi,—X2 ) {5.2b) 

(S"^uj) (xi,X2) = (- 1 )"^UJ ((-l)^Xi,X2 + rmrfn ). (5.2c) 

The glide reflection (or shift-reflect operation) S generates a cyclic group of order 2n, 

C2n = {e,S,Sfl--- 

where e denotes the identity e = S^. The rotation-through- 7 r operation TZ generates a 
cyclic group of order two, R 2 = {e,TZ}. The complete set of discrete symmetries of the 
Kolmogorov equation, therefore, is the dihedral group of order 4n, i.e., 

i74„ = i?2KC'2„ = {e,5,-- - 


Note that the operations TZ and S do not commute, TZS ^ STZ. Instead, we have STZS = 
TZ. 

Therefore, the solutions of the Kolmogorov equation have up to 4n distinct but equiv¬ 
alent copies due to its equivariance under the discrete symmetries Il 4 „. They also have 
infinitely many equivalent copies due to the continuous symmetry 7e for any i £ [ 0 , 27r]. 

An invariant solution may itself have some, all or none of the symmetries of the equa¬ 
tions. The laminar solution Eq for instance has the complete set of symmetries, i.e., 
gEo = Eq for all g £ D^n and Ti.Eq = Eq for all £ £ [0, 27r]. The laminar state, therefore, 
has only one copy. The traveling wave T 7 (see Fig. 8 ), on the other hand, has no sym¬ 
metries and therefore possesses infinitely many equivalent copies. Incidentally, traveling 
wave Tj happens to have the lowest dimensional unstable manifold, dimA“ = 2 , among 
the solutions found here. 

These symmetry related copies greatly complicate the analysis of the state space of 
the Kolmogorov flow. When comparing the if distance between two states and u^, 
one needs to take the minimum distance between and and all their symmetry 
related copies. For example, let u^{t) to be a symmetry copy of a solution u^(t) such 
that u^(t) = (TfU^)(t) for some £ £ (0, 27r). These two states are equivalent and both 
solve the Kolmogorov equation. However, the norm HTfU^ — u^|j 4,2 is generally non-zero. 
Therefore, an appropriate norm on the state space of the Kolmogorov flow is 

min||u^ - , (5.3) 

where the minimum is taken over all £ £ [0, 27r] and g £ D^n- Evaluating this norm can 
be somewhat cumbersome. 
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Another approach is to map each state into a symmetry-invariant polynomial basis. For 
low dimensional dynamical systems with simple discrete symmetries, such coordinates 
are available analytically (Gilmore & Letellier 2007). As the dimension of the system 
(and/or the order of the group) increases, the determination of the invariant coordinates 
becomes quickly prohibitive (Siminos & Cvitanovic 2011). As a result, and to the best 
of our knowledge, symmetry-invariant polynomial coordinates for the Kolmogorov flow 
are not known. 

Here, we take an alternative approach which also proves to be insightful in analysis of 
the temporal intermittency. We define the projection operator 

- 2n—1 

(u + : (5-4) 

An 

m—O 

which is the average over all copies of u given by the discrete symmetries D^n (Cvitanovic 
et al. 2014). The projection T^u is invariant under all discrete symmetries g S and 
therefore we refer to it as the symmetric part of the state u. The symmetric part of 
vorticity w is defined analogously. 

All An symmetry copies of the state u have a unique projection T^u. Working with 
the symmetric part of the states, therefore, eliminates the complications arising from the 
discrete symmetries. 

Considering this symmetric part also has a physical motivation. The energy E, dissi¬ 
pation D and energy input I defined in (4.3) are invariant under symmetry operations. 
Denoting the energy input of a state u by /[u], we have /[u] = /[gu] for g being any 
symmetry of the Kolmogorov flow. The same holds for energy E and dissipation D. The 
particular linear form of the energy input /[u] implies that it is also invariant under the 
projection operation V, i.e. /[u] = I\Pu\. 

Furthermore, we have 

I[u - Vu] = 0, (5.5) 

that is, the remainder u — T^u does not contribute to the energy input. As discussed in 
Section 5.1, the intermittent episodes of the flow are triggered by high energy inputs. 
Therefore, to study the temporal intermittency of the Kolmogorov flow, it is sufficient 
to consider the symmetric part T^u. Note, however, that due to the quadratic form of E 
and D, the energy and dissipation of the remainder u — Vu are generally nonzero. 

The remaining continuous symmetry 7/, is handled here by the method of slices (Car- 
tan 1935; Field 1980; Rowley et al. 2003). This method replaces all continuous symmetry 
copies of a state with a copy that belongs to a given hypersurface called the slice. This 
hypersurface is such that each group orbit 7/u in a neighborhood of a given state inter¬ 
sects the slice transversally at a unique point. Each group orbit 7/u is then replaced by 
its unique intersection with the slice. The method of slices has only recently been used 
in the context of fluid dynamics (Willis et al. 2013, 2015). We use the first-Fourier-mode 
implementation of this method developed by Budanur et al. (2015) (cf. Appendix D for 
more detail). 


5.3. Temporal intermittency 

Figure 12 shows the symmetric part of the vorticity field Vuj for 6 select times along the 
turbulent trajectory. Since the forcing wave number is n = 4, there are 4n = 16 discrete 
symmetries. As a result, the symmetric part of each state exhibits recurring patterns that 
are related through discrete symmetry operations TZ and S. In other words, knowing the 
symmetric part Vui on one-sixteenth of the domain [0, 27r] x [0, 27r], one can reproduce 
Vuj on the entire domain. 



Figure 12. The symmetric part of the vorticity field Vio of the turbulent trajectory at t = 84 
(a), t = 94 (b), t — 104 (c), t = 124.2 (d), t — 134.2 (e) and t = 144.2 (f). The energy 
input/dissipation pairs {I,D) for the states are (0.109,0.107), (0.115,0.114), (0.104,0.104), 
(0.095,0.094), (0.314,0.269) and (0.099,0.093), respectively. 

Panel (e) in Fig. 12, showing the state at time t = 134.2, corresponds to the first 
intermittent peak in Fig. 10(a). A distinct change of topology occurs in the symmetric 
part of the vorticity field as the turbulent trajectory undergoes an intermittent episode. 
Before and after the episode, Vuj has at least two distinct co-rotating vortices in each 
positive (or negative) vorticity band. As the trajectory gets closer to the intermittent 
episode, the co-rotating vortices merge, resulting in horizontal bands of alternating posi¬ 
tive and negative vorticity. After the episode (see panel (f)), these bands split again into 
two distinct co-rotating vortices. This sharp distinction is not immediate from the full 
vorticity field uj (cf. Fig. 13). The same trend (i.e., the merger of co-rotating vortices) 
was observed during the intermittent episodes occurring at later times. 

We found by inspection that the symmetric part of equilibrium solution E 13 (see 
Fig. 14) is strikingly similar to that of the turbulent trajectory as it undergoes inter- 
mittency. This is visually appreciable from comparing Fig. 12(e) with Fig. 14(b). This 
observation suggests that close passages to equilibrium E 13 trigger the intermittent be¬ 
havior. 

In Fig. 14(c), we also show the symmetric part Vuj for equilibrium Eq. While the 
energy input/dissipation of equilibria Eq and if 13 are very close (cf. Table 1), their 
vorticity fields are very different. This demonstrates the fact that closeness of the energy 
input/dissipation of two states does not imply their closeness in the state space. 

To quantitatively examine the role of equilibrium Ei^ on intermittency, we consider 












































































































































Figure 13. The vorticity field ui for the snapshots of the turbulent trajectory shown in Fig. 12. 
The panels correspond to times t = 84 (a), t = 94 (b), t = 104 (c), t = 124.2 (d), t = 134.2 (e) 
and t = 144.2 (f). 
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Figure 14. (a) The vorticity field uj for the equilibrium solution ifia. (b) The symmetric part 
of the vorticity Vlo for the equilibrium solution Fis. (c) The symmetric part of the vorticity Vlo 
for the equilibrium solution Eq. 


the distance between the symmetric part of each turbulent state and the symmetric 
part of equilibrium Ei^, i.e., ||7^a;(t) — 'PEi 3 \\i^ 2 . To account for the continuous transla¬ 
tional symmetry, each symmetric part is first brought to the first-Fourier-mode slice, as 
explained in Appendix D, before the distance is computed. 
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Figure 15. The distance \(Puj{t) — 'Pi ?||^2 between the symmetric part of the turbulent 
trajectory Vuj and the symmetric parts of the equilibria E 13 (top left), En (top right) and 
El (bottom left). The normalized energy input Ijliam of the turbulent trajectory is shown for 
reference (bottom right; same as Fig. 10(a)). Each distinguished peak of the energy input (i.e. 
intermittency) coincides with a significant decrease in the distance from the equilibrium Eiz. 


Figure 15 shows the evolution of the if distance as a function of time. The evolution 
of the energy input / is also shown in the same figure. Every intermittent episode, i.e. 
significant peaks in the energy input, is concurrent with a significant dip in the 
distance. This demonstrates that the intermittent episodes correspond to close passages 
of the turbulent trajectory to equilibrium Ei^. 

The distance between the turbulent trajectory and some of the other invariant solutions 
(e.g. Ell and T 4 ) also decreases during the intermittent episodes. Their decrease is, 
however, not as dramatic as the one observed for equilibrium Ei^. For instance, the 
evolution of the distance from En is also shown in Fig. 15. The distance from En 
fluctuates mostly around 14.0 which is significantly larger that the average distance from 
i?i 3 , which fluctuates around 8.0. During the intermittent episodes, the distance from 
both equilibria decreases. The distance from En reaches 4.64 at its minimum which is 
approximately 1/3 of its average distance of 14.0. The decrease in the distance from 
Ei^ is more significant: it attains a minimum as low as 1.21 which is approximately 6.6 
times smaller than its average distance of 8.0. In contrast, the distance from equilibria 
residing in the ergodic sea increases during each intermittent episode. The distance from 
the equilibrium Ei is shown in Fig. 15 as an example. 

If the Kolmogorov flow at Reynolds number Re = 40 possesses an ergodic attractor, 
as the numerical evidence suggests (see, e.g., Platt et al. (1991)), every invariant solution 
embedded in the attractor will eventually be visited by a generic turbulent trajectory. 
Nonetheless, for the finite-time (and relatively short) turbulent trajectory computed here, 
equilibrium E 13 appears to explain its observed intermittencies. 

Finding the invariant equilibrium and traveling wave solutions is only the first step. A 
complete understanding of the Kolmogorov flow in terms of its invariant solutions will 
require a detailed state space analysis of the type carried out by Gibson et al. (2008) for 
the plane Couette flow. Such a thorough study deserves its own treatment which will be 
presented elsewhere. 













26 M. FARAZMAND 

6. Conclusions and perspectives 

We have proposed and developed here a new method for finding the equilibrium and 
traveling wave (relative equilibrium) solutions of the forced Navier-Stokes equations with 
periodic boundary conditions. Namely, adjoint partial differential equations (PDEs) were 
derived whose (relative) equilibria include those of Navier-Stokes equations. Furthermore, 
the (relative) equilibria of the adjoint equations are asymptotically stable, and therefore, 
their trajectories converge to desired invariant solutions. 

Applying this method to the Kolmogorov flow led to the discovery of several new 
equilibrium and traveling wave solutions. Specifically, for Reynolds number Re = 40, we 
found a total of 24 non-trivial equilibrium and traveling wave solutions, where only 3 of 
them were known previously (Chandler & Kerswell 2013). 

Some of these new solutions exhibit highly localized spatial structures (cf. figures 6 
and 8). Such localized structures were previously believed to only exist in rectangular 
domains with large aspect ratios (Schneider et al. 2010; Lucas & Kerswell 2014). 

One of the equilibrium solutions appears responsible for the observed temporal inter- 
mittencies in the Kolmogorov flow. Such intermittencies manifest themselves as sudden 
and short-lived bursts in the energy dissipation (and energy input) along a generic turbu¬ 
lent trajectory. We showed that these bursts correspond to close passages of the trajectory 
to equilibrium ifia. This supports the dynamical systems view that such rare, extreme 
events are the result of heteroclinic orbits carrying the trajectory to less accessible regions 
of the state space (Holmes 1993). 

Due to the periodic boundary conditions, the adjoint PDEs can be easily integrated 
numerically using a pseudo-spectral method. The spectral representation of the adjoint 
equations, presented in Appendix C, exhibits a close resemblance to that of the Navier- 
Stokes equations. Therefore, existing pseudo-spectral codes can be easily adapted to solve 
the adjoint equations. 

We found the rate of convergence of the adjoint method to be rather slow. More 
precisely, the trajectories of the adjoint PDEs converge rapidly to the vicinity of an 
invariant solution, but further convergence takes place at a slow rate. To achieve better 
convergence, we proposed a hybrid adjoint-Newton algorithm, consisting of two steps: 
First, the adjoint equations are integrated from some initial condition in order to reach 
the neighborhood of an invariant solution. Once in the neighborhood of the invariant 
solution, the standard Newton-GMRES-hook iterations (Viswanath 2007) were used to 
converge further to the solution. This hybrid algorithm yielded a 100% converge rate 
from generic initial conditions (4.6), obviating the preprocessing step for finding ‘good’ 
initial guesses (Viswanath 2007; Chandler & Kerswell 2013). 

Newton-GMRES-hook iterations are relatively expensive. Hence, a modification of our 
adjoint-based method, that would eliminate the Newton-type steps altogether, is highly 
desirable (see Lakoba & Yang (2007), in the context of nonlinear wave equations). 

While we only considered periodic boundary conditions, the general adjoint-based ap¬ 
proach applies to wall-bounded turbulence such as plane Couette and pipe flows. Our 
preliminary analysis (not presented here) shows that, in the presence of boundaries, the 
resulting adjoint equations require more boundary conditions than the corresponding 
Navier-Stokes equations. This is to be expected as the adjoint equations have higher 
order spatial derivatives. This calls for a special numerical treatment of the adjoint equa¬ 
tions for wall-bounded flows. Our results on the Kolmogorov flow, however, show that 
the gain is worth the pain. 

Finally, we point out that our adjoint-based method does not immediately apply to the 
computation of periodic and relative periodic orbits. They are currently found through 
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Newton-GMRES-hook iterations (see, e.g., Kawahara et al. (2006); Viswanath (2007); 
Willis et al. (2013); Lucas & Kerswell (2015)) with the drawbacks discussed in the Intro¬ 
duction. Alternatives include the variational method of Lan & Cvitanovic (2004) which 
shares the universally convergent property of our adjoint-based method. Its computa¬ 
tional complexity is, however, formidable (Fazendeiro et al. 2010). More recent method 
of Yang (2015) has proved promising for unidirectional wave equations, but its feasibility 
for Navier-Stokes equations remains to be explored. 
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A tutorial, with the accompanying MATLAB code, is available at 
https: / /farazmand. wordpress. com/2015/12/18/ adj oint / . 


Appendix A. Derivation of the adjoint equations for equilibria 

We start by redefining the vector Fq by adding the term V • u as a new element; that 
is 

+ (Al) 

where q = (u,p) and with the understanding that for divergence-free vector fields, the 
last row of Fq is identically zero. This twist in the notation proves useful below, where 
we compute the adjoint of the Gateaux derivative within the space of divergence-free 
vector fields. The Gateaux derivative of Fq is given by 


^o(q;q') 


—u' • Vu — u • Vu' — Vp' + z/Au'\ 

V -u' j 


(A 2) 


where q' = (u',p') with V • u' = 0. 

The key part of deriving the adjoint equation (3.3) is to find the adjoint operator jCq. 
Its derivation is standard and can be found in the literature on adjoint-based optimal 
control (see, e.g., Gunzburger (2002); our notation is closer to Farazmand et al. (2011)). 
The difference is that, in optimal control, one seeks to minimize a cost functional with 
the constraint that the Navier-Stokes equations are satisfied. Here instead, the only 
constraint is the divergence-free condition and we seek to find the states q = (u, p) such 
that ||Fo(q )||.4 = 0. We derive the adjoint with respect to the inner product first. 
The adjoint with respect to the more general inner product (2.2) follows easily. 

Let the function space H be the space of square integrable functions q = (u,p) such 
that the u component is divergence-free. More precisely, 

•H = {q = (u,p) e L'^iV) I V • u = 0}. 


(A3) 
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We define the usual inner product on H, i.e., for any q, q' S H we define 

(q,q')L 2 =/' [u(x,t) • u'(x,t)+p(x,ty(x,t)]dx (A4) 

Jv 

The adjoint of the linear operator (A 2) with respect to the inner product (A 4) on the 
function space H (A 3) is given by 


^o(q;q") = 


[Vu" + (Vu")^] U - Vp" + 

V-u" 


Vq, q" e H, 


(A 5) 


where T denotes matrix transposition. As we are restricted to the space of divergence free 
vector fields, the last component of the adjoint (A 5) is identically zero, i.e., V • u" = 0. 

This expression follows directly from definition (2.11) of the adjoint and Eq. (A 2). 
Namely, substituting (A 2) into (2.11), we have 


(£o(q;q'),q")L= 


/V 


—u' • Vu — u • Vu' — + iyAu'\ ( u"\ , 

' ' » ) dx. 


V • u' 


V 


= / [((Vu" + Vu"T) u - Vp" + z/Au") • u' + (V • u")p'] dx, 
JT> 

(Vu" + Vu"^) u - Vp" + i^Au"\ fu‘ 

V-u" ) ■ i p' 


/V 


dx. 


= (q'yo(q;q")) 




(A 6) 


where the first line is the definition of the inner product and the second line follows 
from integration by parts. Note that the boundary terms vanish due to the periodic 
boundary conditions. Only one of these integration by parts is not straightforward, which 
we outline below. 


(-u' • Vu 


/X) 


u • Vu') • u" dx = 


L 

L 

L 

L 


{—u”u'jdjUi — u'lujdju'j) dx, 
{uidj(u'-Uj) + u[dj{u”uj)) dx, 
{uiu'jdju'l + u'^Ujdju'l) dx, 
{uju[diu” + u[ujdju”) dx. 


= / [(Vu" 


Vu"T) ul 


u' dx, 


/x> 


where we used the summation notation on repeated indices and the fact that u and u' 
are divergence-free: djUj = dju'j = 0. 

Identity (A6) holds for general q, q',q" € T-L. Along the trajectories of the adjoint 
descent PDE, we have q = (u,p), q' = {drU,drp) and q" := (u",p") = Fo(q) (cf. 
Eq. (2.10)). This yields. 


(q',£|,(q; q"))^^ = / {^u • ([Vu" + (Vu")^] u - Vp" + ^Au") + (V • u")yp} dx. 
Jv 

(A7) 

Therefore, equation (3.3a) ensures that the above inner product is negative semi-definite. 
This together with identity (2.10) ensures that the norm |jFo(q(T))||i 2 decreases along 
the trajectories of the adjoint PDE. 

Due to the divergence-free condition, the term (V • u")dTP vanishes. Therefore, no 
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independent evolution equation is obtained for the pressure p. As in the case of the 
Navier-Stokes equation, the pressure terms simply enforce the divergence-free condi¬ 
tions (3.36). 

The adjoint descent with respect to the more general inner product (2.2) is obtained 
similarly, the only difference being that q" is replaced by Acf'. 


Appendix B. Derivation of the adjoint equation for traveling waves 

As in Appendix A, we redefine Fc as 

—u • Vu — Vp -1- i/Au -I- f -I- c • Vu 

Vu 

where q = (u,p) and 

Fc(q)=Fo(q)+(^‘'’^"), (B 2) 

with Fq defined in (A 1). Due to the linearity of the second term, the Gateaux derivative 
of Fc is 

^c(q;q') =>Co(q;q') + 

where q' = (u',p') and Co is defined by (A 5). 

Analogous to anlysis of Appendix A, the adjoint can be shown to be given by 

^c(q;q") ='Co(q;q") - 

where q" = (u",p"). 

We would like to define the adjoint descent equation in such a way that, along its trajec¬ 
tories (q(T), c(t)), the residue ||Fc(i-) (q(''‘))IU decreases monotonically. Taking derivative 
with respect to the hctitious time t, we obtain 

d.||Fc(q)||5, = (i:c(q;q'),^Fc(q))L2 + ^(^^-^"^ ,AFc(q)^^^ 

= (q',>Cl,(q; AFc(q)))i2 -t- (c • Vu,u")^2 

= (q',>Cl,(q;AFc(q)))i 2 -He • /" (Vu)Tu"dx, (B5) 

Jv 

where q' = 9rq and u" is given by (3.10). The adjoint set of equations (3.9) is designed 
in such a way that the right hand side of the above equation is always negative, resulting 
in monotonic decrease in the residue along its solutions (u(T),p(r)). 




Appendix C. The adjoint equations in the Fourier space 

In this Appendix, we derive the spectral representation of the adjoint descent equa¬ 
tions (3.3) and (3.9). As shorthand notation, we define 

N := [Vh"-H (Vh")^] u, (Cl) 

and denote Fourier transforms with a hat sign. Then one can write equation (3.3a), in 
the Fourier space as 
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where I is the identity matrix and k (g) k denotes the tensor product. Note that the 
pressure p" is eliminated using identity (3.6) which implies 


p"(k) = 

For u", using definition (3.11), we have 


-ik-N(k) 

■ 


u"(k)= 


l + |k|2- 

The term u" is in turn computed from identity (3.4) as 

kgk 


u"(k) = - I- 


|k|^ 


R(k) - v\k.f u(k) + f(k), 


where R is the shorthand notation for the nonlinear term 

R := u Vu. 


(C3) 


The nonlinear terms R and N are computed by the pseudo-spectral method, i.e., the 
differentiations are carried out in the Fourier space and the products are computed in 
the physical space. 

Note that the spectral representation (C 2) of the adjoint equation closely resembles 
that of the Navier-Stokes equation (C3). Therefore, existing pseudo-spectral codes for 
Navier-Stokes can easily be adapted to solve the adjoint PDFs 

The adjoint descent for the traveling waves is computed similarly, except that the term 
i (c • k)u"(k) is added to the right-hand-side of equation (C2) accounting for the term 
c • Vii" in (3.9). More precisely, the adjoint descent equation for the traveling waves in 
the Fourier space reads 

^ N(k) -H z/|kp fi"(k) -h i (c • k)a"(k), (C4) 

with u" = An" and 

R(k) - z/|kp u(k) -I- f(k) -I- i (c • k)u(k). (C 5) 


u"(k) = - I- 


|k| 



Appendix D. Symmetry actions in the Fonrier space 

The symmetry operations (5.1) can be readily implemented in the Fourier space. The 
following transformations follow directly from the definition of the Fourier transform. 
Denoting the Fourier modes of the velocity field u by u(k) with k = (fci, ^ 2 ) G we 
have 

u(xi,X2) = 51 H u(fci,fc2)e*'=^"^e*'=="T (D 1) 

felGZ fc2GZ 

Therefore, the shift operation Ti satisfies 

{Tin){xuX2)=n{x^+i,X2) = ^ E 2(^1, ^^ 2 ) 6 *'=^ (D2) 

fciGZ fe2GZ 

implying that the action of the continuous symmetry 7^ on the Fourier mode u(fci,fc 2 ) 
is multiplicative, such that 


7^u(fci, fca) = ^ 2 ). 


(D3) 
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Similarly, for the rotation through tt symmetry TZ, we have 

7^u(fcl, ^ 2 ) = -u{-ki,-k2) = - [u{ki,k2)]* , (D 4) 


where the superscript * denotes complex conjugation and the last identity follows from 
the fact that u is real valued. 

Finally, the shift-reflect symmetry S acts on the Fourier modes through 


S^u{ki,k2) = 


f {-irui{{-l)^ki,k2)\ 

V M2 ((-l)’"fcr,^2)/ ■ 


(D5) 


To reduce the continuous symmetry, we use the first-Fourier-mode slice of Budanur 
et al. (2015). It follows from (D 3) that the shift Te acts on the Fourier mode (^ 1 ,^ 2 ) = 
(1,0) of the vorticity through 

w(l,0) w(l,0)e*^ = |w(l,0)|e*('^(^’°)+^^ (D6) 


where 4>{ki,k2) G (—7r,7r] denotes the principal value of the phase of mode (^ 1 ,^ 2 ). 
To bring the vorticity field to the hrst-Fourier-mode slice, the shift value i is chosen 
such that ^i(l,0) + £ = 0. More precisely, a given vorticity field uj is replaced with its 
symmetry-equivalent copy through the transformation 

Ujiki,k2) ^ A:2). (D 7) 


As a result, the mode w(l, 0) of the symmetry-reduced vorticity field has vanishing imag¬ 
inary part. 
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